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Abstract 

This paper deals with stability in the numerical solution of the promi- 
nent Heston partial differential equation from mathematical finance. We 
study the well-known central second-order finite difference discretization, 
which leads to large semi-discrete systems with non- normal matrices A. 
By employing the logarithmic spectral norm we prove practical, rigorous 
stability bounds. Our theoretical stability results are illustrated by ample 
numerical experiments. 
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1 Introduction 



This paper deals with stability in the numerical solution of the Heston partial 
differential equation (PDE), 



du -, 2 d 2 u d 2 u i 2 d 2 u du du 

m =2SV d^ + P(TSV d^ +2<JV W + "Ta + «* - v) \Tv 



for s > L, v > and < t < T. The Heston PDE constitutes one of the 
prominent equations of mathematical finance, cf. e.g. [H HI [9j [11]. It gener- 
alizes the celebrated one-dimensional Black-Scholes PDE where the volatility 
is modelled by a stochastic process rather than being constant. Clearly, (jl.ll) 
can be viewed as a time-dependent advection-diffusion-reaction equation on an 
unbounded two-dimensional spatial domain. The exact solution value u(s,v,t) 
represents the fair price of a European-style option if at time T — t the under- 
lying asset price and its variance equal s and v, respectively, where T > is 
the given maturity time of the option. The quantity L > is a lower barrier, 
K > is the mean-reversion rate, rj > is the long-term mean, a > is the 
volatility-of- variance, p £ [—1,1] is the correlation between the two underlying 
Brownian motions, and r > is the interest rate. These quantities are all given 
and arbitrary. We remark that in practice the correlation p is usually nonzero, 
and hence, (|1.1[) contains a mixed spatial-derivative term. The Heston PDE is 
complemented with initial and boundary conditions which are determined by 
the specific option under consideration. In this paper we shall assume boundary 
conditions of Dirichlet type. 

A widely known semi-discretization of PDEs in finance is given by central 
second-order finite difference (FD) schemes, see e.g. |11[ IT3] . To render the 
numerical solution of the Heston PDE feasible, the spatial domain is first re- 
stricted to a bounded set [L, S] x [0, V] with fixed values S, V chosen sufficiently 
large, with additional Dirichlet conditions imposed at s = S and v — V. Let 
mi , mi > 3 be any given integers and define spatial mesh widths 

S ~ L a V 
As = , Av = 



mi + 1 m 2 + 1 

The central second-order FD schemes for approximating the advection, diffusion 
and mixed derivative terms in (|l.lj) arc 

(u v ) u « ''' J (L2b) 



2Av 

Ui+ij - 2ujj + Ui-ij 

jAsf 
Uij+i - 2ujj + Uij-i 

Ja^Y 2 



(u ) w "" l ' J+1 - ^ u,t,J ~ 1 (1 2d) 



4AsAu 



with the short-hand notation Uij — u(si,Vj,t) and spatial grid points 

Si — L + i ■ As (i = 0, 1, . . . , mi + 1) , Vj = j ■ Av (j = 0, 1, . . . ,ma + 1). 

Semi-discretization by (jl.2[) of a given initial-boundary value problem for the 
Heston PDE leads to an initial value problem for a large system of ordinary 
differential equations (ODEs), 

U'{t) = AU(t) + b(t) (0<t< T), U{0) = U Q . (1.3) 

Here A is a given constant real m x m matrix and b{t) (for < t < T) and Uq 
are given real m x 1 vectors with m = m\mi. The vector Uq is directly obtained 
from the initial condition for whereas the vector function b depends on 

the boundary conditions. For each t > 0, the entries of the solution vector 
U(t) to (|1.3[) form approximations to the exact solution values u(si,Vj,t) for 
1 < i < mi, 1 < J < ni2. 

The aim of our paper is to gain insight into the stability of the semi-discrete 
Heston PDE (|1.3p . To this purpose, we are interested in the existence of useful, 
rigorous upper bounds on the quantity ||e <y!l || (for t > 0) where || • || denotes an 
induced matrix norm. Such bounds, on the magnitude of the matrix exponen- 
tial of tA, guarantee that any (rounding or discretization) errors cannot grow 
excessively. For central second-order FD discretizations of the Black-Scholes 
PDE, adequate stability bounds were recently proved in [6]. These bounds are 
of the well-known type 

\\e tA \\ < Ke tu (t>0) (1.4) 

with constants w € K and K > 1. To our knowledge stability estimates of the 
type (jl.4l) have not been obtained in the literature up to now for FD discretiza- 
tions of the Heston PDE. In the present paper, we shall establish a natural 
extension of stability results derived in [B]. We note that a main difficulty in 
proving this extension lies in the mixed derivative term in the Heston PDE, 
which does not arise in the Black-Scholes case. 

As the semi-discrete Heston matrix A is in general non-normal, bounds on 
the norm of e tA which are based solely on the eigenvalues of A are most often not 
useful. For the stability analysis in this paper, we shall employ the logarithmic 
spectral norm. For any given complex k x k matrix A, with integer k > 1, it is 
defined by the limit 

||7 + tA|| 2 -l 
(i 2 [A\ = hm , 

tio t 

where || • 1 1 2 is the spectral norm and I is the k x k identity matrix. We note 
that general complex matrices A are considered for later use. The following key 
result forms the basis for our analysis; see e.g. pi \7\ \TU \ fT2 ] . 

Theorem 1.1 Let A be any complex k x k matrix and uj £ R. Then 
(i 2 [A}<co \\e tA \\ 2 < e tuJ for all t > 0. 
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Denote by (• , -)2 and | • I2 the standard inner product and Euclidean norm, re- 
spectively. Then for the logarithmic spectral norm one has the more convenient 
formulas 

^ 2 [A] = max{Re(Ax,x) 2 : x £ C k , \x\ 2 = 1} (1.5a) 
= max {A : A eigenvalue of \{A + A*)} , (1.5b) 

where A* stands for the Hermitian adjoint of A. 

Motivated by the study [6] for the Black-Scholes PDE, we introduce also a 
suitably scaled version of the spectral norm on C rrixm . Consider the positive 
diagonal matrices 

Di = diag(si,s 2 ,...,s mi ), D 2 = diag(ui,u 2 , . . . ,v m2 ) , D = D 2 ® D x , 
where <g> is the Kronecker product. For vectors x £ C m we define the norm 

\x\ D = \D-^ 2 x\ 2 

and denote for matrices A £ (£™xm. ^ e mc j uce( j matrix norm and logarithmic 
norm by ||^4||_d and /z_d[A], respectively. For any matrix A there holds 

\\A\\ D = \\D- 1 ^AD 1 / 2 \\ 2 , ^ D [A]=n 2 [D- 1 l 2 AD 1 ' 2 ] (1.6) 

and the spectral norm of A is bounded in terms of its scaled version through 

V s i v \ 

The outline of the paper is as follows. In Section [2] we derive practical 
stability bounds for the semi-discrete Heston PDE (|1.3I) . Here the advection 
and diffusion terms are each studied individually. Numerical illustrations are 
provided in Section [3l with actual computations of the norms of matrix expo- 
nentials. Conclusions and issues for future research are discussed in Section |U 

2 Stability bounds 

Let / denote the identity matrix of generic dimension. Associated with the FD 
formulas (|1.2j) . we define the tridiagonal m\ x mi matrices 

Li = 2^-tridiag(-l,0,l) , M x = ■ tridiag (1 , -2 , 1) 

and the tridiagonal m, 2 x m 2 matrices 

U = ^j~- tridiag (- 1, 0, 1) , M 2 = j—^ ■ tridiag (1 , -2 , 1) . 
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FD discretization by (|1.2[) of the spatial derivative terms rsu s , k(jj — v)u v , 
^s 2 vu ss , pasvu sv , ^a 2 vu vv in the Heston PDE (jl.ip gives rise to the following 
real m x m matrices, respectively: 



Ax = r/ (2.1a) 

A 2 = K[{r}I-D 2 )L 2 ]®I, (2.1b) 

A 3 = \D 2 ® , (2.1c) 

At = po{D 2 L 2 ) ® (DiLi) , (2.1d) 

A 5 = \a 2 {D 2 M 2 )®I . (2.1e) 



Here a lexicographic ordering of the spatial grid points is considered. It is worth 
noting that (D 2 L 2 )(g> (DiLi) in (|2.1H ) can be regarded as a discrete analogue of 
(to„) o (su s ) = svu sv where o denotes composition. The semi-discrete Heston 
matrix A in (|1.3[) is equal to 

A = Ai + A 2 + A 3 + A 4 + A 5 - rl. 

Our introductory result concerns the two parts of the semi-discrete Heston 
matrix corresponding to the advection terms in the s- and w-directions. It 
provides useful stability bounds of the type (|1 .4[) for these. 

Theorem 2.1 Let r,K,rj > and let A\, A 2 be given by \2. lb ). 112. lb ). Then 
\\e tAl \\ 2 < e tu (t > 0) with u= T - 

and 

\\e tA2 \\ 2 < e tuJ (t>0) with u = -. 

The above values for uj are the smallest that hold uniformly in the respective 
mesh widths. 

Proof Consider the symmetric matrix F = tridiag (i, 0, i) of generic dimen- 
sion. All eigenvalues of this matrix lie in the real interval [—1,1]. It is readily 
verified that A\ + Aj = -rl ® F and A 2 + A^ = kF ® I (for any 77). The 
eigenvalues of I ® F and F ® I are the same as those of the pertinent matrix F. 
By p. 5b ) it thus follows that /^[Ai] < r/2 and /Lt2[A2] < k/2 and application of 
Theorem 11.11 yields the required bounds. Furthermore, as there exist eigenval- 
ues of F that converge to —1 and 1 when the dimension increases, the obtained 
values for ui are the smallest that hold uniformly in the respective mesh widths. 

□ 

The subsequent lemma deals with the logarithmic spectral norm of certain 
matrices of block Toeplitz type and is essential to the proof of our main result 
in this paper. Let E = tridiag(0, 0, 1) denote the m 2 x m 2 forward shift matrix. 



5 



Lemma 2.2 Let Bq, B\ be any given real mi x m\ matrices and let the to x to 
matrix B be defined by 

B = I®B +E®B 1 +E T ®Bj. 

Then 

fx 2 [B]< max fi 2 [B + 2C,B X ]. 
CeC,|C|=i 

Proof Consider the so-called symbol of B, given by 

B(Q =B + (B 1 +C 1 Bj 
for ( e C, |C| = 1. Since B is a block Toeplitz matrix, also the exponential 

P 



e' B 

j=0 ■' 



is block Toeplitz. The symbol of e tB is equal to e tB ^ and one has the bound 

||e tB || 2 < max ||e tB ( c >|| 2 , 

C€C,|C|=1 

which is a consequence of Parseval's identity, see e.g. [T[ p. 186]. By Theorem 
II. H it readily follows from this that 

MB] < m«a M£(C)]. 

Let i?(C) = Bq + 2(Bi. Then the Hermitian parts of B(() and B(() are equal 
and hence, by (|1.5b ). there holds fi 2 [B(C)] = \i 2 [B{()}- This yields the proof. 

□ 

Our main result of this paper concerns the stability of the diffusion part 
(including the mixed derivative term) of the semi-discrete Heston system. 

Theorem 2.3 Let a > and p £ [—1, 1] and let A3, A4, A$ be given by A2.1k ), 
fOH). fOk). Then, for all t > 0, 

\\e t{A3+M+M) \\ D < 1, (2.2a) 

J(A 3 +A 4 +A a )r IS m ,Vr, 



n e t( A3+Ai+A5) n < m . (2.2b) 
V 

The strong stability result (|2.2a ) means that the diffusion part of the semi- 
discrete Heston system is contractive in the scaled spectral norm. The bound 
(12.2b ) for the standard spectral norm is discussed in more detail in Section [3J 
Theorem l2.3l can be viewed as a natural extension of [6j Theorem 2.8] that was 
derived for the case of the Black-Scholes PDE. In the special situation where 
p — 0, so that no mixed derivative term is present in the Heston PDE and the 



G 



matrix A4 vanishes, the result of Theorem [531 can be obtained in analogous way 
to loc. cit. However, the important general situation where p ^ requires a 
new, and more elaborate, proof. 

Proof The bound pb) follows directly from ([2~2"h ) by (TO]). By Theorem [TTTT 
the bound (|2.2fa ) is equivalent to /z_dL4 3 + A 4 + A 5 ] < 0. In the following we 
show that this condition holds. For convenience, the proof is split into three, 
consecutive parts. 



(i) For any given real square matrices A, G with G nonsingular it holds that 
fi 2 [A] < if and only if p 2 [G T AG] < 0. Choosing A = A 3 + A A + A 5 and 



G = D. 



-1/2, 



, we obtain 
=> ^[B}<0, 



/, and taking into account 
p D [A 3 + A 4 + A 5 ] <0 
where the matrix B is given by 

B = (D2 1 ®D^ 1/2 )(A 3 +A 4 + A 5 )(I®D 1 1 /2 



\1 <g> {D\l 2 MiD[ /z ) + poL 2 ® {D^LxD 1 ^) + \a 2 M 2 <g> I 



,1/2 



1/2, 



l 1 ^ 



Let a = a/Av and define the matrices 



Li = L^ /2 Li£>{ /2 , Mi = D\ /2 MiD\ 12 . 



Note that 



La 



1 



2Av 

Inserting into B yields 



(E £ T ) 



M 2 



(Ay)' 



T (E-2I + E T ). 



B = 5I ® Mi + ® Li + 5(T 2 M 2 <8> J 



J® Mi + pa(E - E T ) ® Li + 65 2 (£ - 21 + E T ) (g> I 
I® (Mi- 2d 2 1) + E® (paLx + a 2 1) + E T ® (-paLi + a 2 1) 



Since Lj — —L\ we are in the situation of Lemma 12.21 Application of this 
lemma yields the following sufficient condition for n 2 [B] < 0, with £ £ C : 



//■2 



±Mi - a 2 1 + ((paLi + a 2 1) < whenever |C| = 1. 



Let i denote the imaginary unit and let A max [A] stand for the maximum eigen- 
value of any matrix A having just real eigenvalues. Using (11.5b ) one readily 
finds that the sufficient condition above is equivalent to 



An 



i(Mi + Mi) + 2i(ImC)|o5Li < 25 2 (l-ReC) whenever |C| = 1- (2.3) 
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(ii) Define C s = D\L\ and C ss — \D\M\. Remark that these matrices can be 
viewed as FD discretizations of the su s and ^s 2 u ss terms, respectively. Clearly, 

Lx = D~ 1/2 C S D\ /2 . 
Next, a direct calculation shows that 

|(M 1 £>i-£>iM 1 ) = Li 

and using this one obtains 

|(Mi + mT) = d; 1/2 (2c ss + c s )d\ /2 . 

Therefore, by a similarity transformation, (|2.3p is equivalent to 

A max [C ss + \C S + i(ImQpaC s ] < a 2 {\ - RcC) whenever |(| = 1. (2.4) 
For |C| = 1 it holds that 

1 -RcC> |(l + ReC)(l-ReC) = \{l - (ReC) 2 ) = ±(ImC) 2 . 
This bound gives rise to the following sufficient condition: 

A max [Css + \C S + iypvC s ] < \a 2 y 2 whenever y G M, \y\ < 1. 

Then, upon replacing ^y pa by y and using that the correlation p satisfies \p\ < 1, 
we arrive at the neat condition 

A max [Css + (| + 2iy) C,] < 2y 2 whenever y e R. (2.5) 

Summarizing, 

(E3D =>> (El (EH) =► ^[£]<o. 

In the third and final part we prove that (|2.5[) is fulfilled. 

(iii) Let /iooL4] denote the logarithmic maximum norm of any complex square 
matrix A. It is well-known that if A = (a^j) then 

Poo [A] = max ( Re a*; + |a; 1 ) . 

i z — ' 

Any induced logarithmic norm forms an upper bound on the real parts of the 
eigenvalues of A. In the following, the logarithmic maximum norm will be used 
to this purpose. 

Write Vi — S{/As for 1 < i < m%. There holds 

C ss + (f + 2iy) C s = tridiag(ft, a u 7,) 
8 



with 



oti = -v} , /3 t = \vi (yi-\- 2iy) , 7 l = \v t (i>i + | + 2iy) . 

In proving (|2.5p we need to distinguish two cases: \y\ > 1/2, and the more 
intricate case |y| < 1/2. 



1 2/| > 1/2 Put Pi = 0, 7„ H = 0. One has 



A max [C ss + (| + 2iy) C,] < Moo [C ss + (I + 2iy) C s ] = ^ max ^ai+lftl+^l}. 
Let 1 < i < mi. Then a t + + |-y € ] < 2y 2 if 



where 9 = 4y 2 . By an elementary calculation one verifies that this inequality is 
equivalent to 

40(0 - l)vf + 9 2 (46 - + 9 4 > 0, 
which holds whenever > 1. Thus, condition (|2 .5[) is valid whenever |y| > 1/2. 



|jy| < 1/2 Let A = diag(<5i, (52, • • ■ ,S mi ) with arbitrary real numbers 8{ > 



(1 < i < mi) and write £j = Si/5i-i (2 < i < mi). A similarity transformation 
with the diagonal matrix A leads to the following bound, 

A max [C ss + (| + 2iy) C B ] < Moo [A (C ss + (± + 2iy) C s ) A" 1 ] = 

maxiai + — |7i| , max {a<j+£i|/%|H — 1 7 -|} ; a + e \p 

L £2 2<i<nn-l e i+ \ 



Let 2 < z < mi — 1. The estimate 



yields 



where 



2y 2 



\x±2iy\<x+ — (for x G E, x > 0) 



o» + £i\/3i\ + — !— |7i| < a l + 6j ■ 2y 2 , 



ft, 



6* = 



2 
2 



-2l/< + e<(l/« - |) + + I) 

£i+l 



Ei+l(i/j + 2) 



It holds that 

a; + bi ■ 2y 2 < 2y 2 (whenever \y\ < |* 



aj < and 2a,i + bi < 1 . 
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By trial and error, we have found the convenient choice 



(for 2 < j < mi). 



(2.6) 



Using (I2.6p . and noticing that v i+ i =1^ + 1, it is easily seen that 



1 v, 



a 1 



« 4 



8^0, + |) 



Since Vi>i>2 there follows cij < 0. Next, 

i , (^ + 1) 2 



v- 



[Vi + ±) 2 (^ + §) 



A straightforward calculation shows that 



,3 3, ,2 3, 



-2- > 

16 — u - 



It is readily verified that the inequality in the right-hand side is fulfilled. Hence, 
with flUj) , 

a* + Ei|AI + |7i| < 2y 2 (2<i<mi-l). 

By an analogous reasoning it follows that 

on + — |7i| < 2y 2 , a mi +e TOl |/3 mi | < 2y 2 . 
£2 

Consequently, condition (|2.5p is also valid whenever \y\ < 1/2. This completes 
the proof of the theorem. 



□ 



3 Numerical experiments 

In this section we numerically examine the stability bound (|2.2b ) of Theorem 
I2.3l for the diffusion part of the semi-discrete Heston system. Slightly rewritten, 
it reads 

|| e t(A,+iU + A.)|| < + 

V m\L + b 

The right-hand side is equal to \Jnv\m% = y/m if L = 0, and it is at most equal 
to y/mm{mi,S/L} ■ mi whenever L > 0. For L > the stability bound (|2.2b ) 
is thus more favorable than for L = 0. 

We estimated in MATLAB (version R2009a) the maximum of \\e t(A3+Ai + A ^\\ 2 
over t > for a variety of cases. We considered all combinations of parameter 
values 

a £{0.1,0.2} , pe {-1,0,1} , Le{0,10}. 
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Figure 1: Graph of estimated max t > ||e*( A3+A4+A ^||2 vs. m 2 = 5, 7, 9, . . . , 25 
for L = (black squares) and L = 10 (grey circles) where mi = 2m 2 - Left 
column: er = 0.1. Right column: a = 0.2. Top row: p = 1. Middle row: p = 0. 
Bottom row: p = — 1. 
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Following [5] we chose mi = 2m,2 (so that the dimension m = 2777,2) and se- 
lected m 2 = 5, 7, 9, ... , 25. Further S = 800, V — 5 were taken as in loc. cit. 
For the computation of the matrix exponential and the spectral norm the MAT- 
LAB functions expm and norm (-,2) were used. We note that the feasibility of 
expm implied ni2 = 25 as the largest reasonable choice (then m = 1250). The 
maximum over t > was estimated in a basic way by sampling the values for 
t — 0, 1,2,. ..,100 and subsequently refining in the region around the largest 
value. We mention that the location of the maximum was always found to lie 
in the interval < t < 5. 

The obtained results are displayed in Fig. [TJ Each of the six subfigures 
shows the estimated maximum of ||e t (" 43+j44+ ' 45 )||2 over t > versus 7712 for a 
given pair (a, p) . The black squares correspond to L = and the grey circles to 
L — 10. As a first observation from Fig. [IJ it is readily seen that all numerical 
results are in agreement with the theoretical stability bound (12.2b ). Secondly, 
Fig. [1] reveals that for L — 10 the computed maximum of ||e^ j43+A4+A5 '||2 is 
never larger, and in general much smaller, than that for L = Q. In addition, 
we find in all cases a growth that appears to be at most directly proportional 
to y/m ~ m-2 and ^/to2 when L = and L — 10, respectively. This agrees 
with the bound (|2.2b ) as well, as discussed above. Thirdly, Fig. Q] indicates the 
positive result that the value of a and especially p only has a limited impact on 
the actual maximum of | |e*^ 3+ ^ 4+ ^ 5 ' | | 2 . Note that for p we considered here the 
interesting extreme cases —1, 0, 1, but this result was confirmed by numerical 
experiments with various other values. 

4 Conclusions and future research 

In this paper useful, rigorous stability bounds have been derived relevant to 
central second-order finite difference discretizations of the Heston PDE from 
mathematical finance. Results for the advection and diffusion parts have been 
proved individually and are valid for arbitrary Heston parameters. The stability 
estimates obtained in this paper can be viewed as natural extensions of recent 
stability results from [5] for the case of the one-dimensional Black-Scholes PDE. 

Besides the standard spectral norm, a suitably scaled version has been con- 
sidered, following a fruitful idea from loc. cit. The main result of our paper, 
Theorem [231 states that in this scaled spectral norm the semi-discrete diffusion 
part of the Heston PDE is contractive. This result holds for arbitrary correla- 
tion values p e [—1,1] and thus covers the practically important situation where 
a mixed spatial-derivative term is present. 

The bound in the standard spectral norm is (also) uniform in p, which has 
been illustrated by ample numerical experiments. Both theoretical and numer- 
ical evidence reveals that in the standard spectral norm the stability of the 
semi-discrete diffusion part is much more favorable if the lower barrier L > 
than if L = 0. In actual applications, L > is often fulfilled, for example for 
barrier options; else it is harmless to increase L slightly, when the actual region 
of interest for the asset prices lies far away from this value. 
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We note that the results in this paper can directly be combined, using a 
well-known theorem due to von Neumann [3l Sects. IV. 11, V.7], to arrive at 
stability bounds for various classes of time-discretization schemes applied to 
the semi-discrete Heston PDE, e.g. Runge-Kutta methods and linear multistep 
methods. For the sake of brevity we have not explicitly included these results 
here. 

In future research we shall investigate, among others, the stability of FD 
schemes for the Heston PDE on non-uniform spatial grids. Such grids play an 
important role in mathematical finance. In [5J 113) stability bounds pertinent 
to non-uniform grids were derived for the case of the Black-Scholes PDE and 
more general one-dimensional advection-diffusion-reaction equations. In future 
research we also intend to study for example the adaptation of the obtained 
stability results to different types of boundary conditions. 
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